%%%%---- FIGURE D.1 -----%%%%

%% Plot DELTA and costs
clc;
clear all;
%%
% use section4_all, has 100% of observations
load('section4_all.mat');
load('lambdas_hat.mat');
format short g
%% 
% unpack data needed
donation=d.donation; default=d.default;
relfreq=d.relfreq;
lambdas=x';
n=length(relfreq);
% calculate Delta: gain from opting out
Delta=Delta_fun(donation, default);

%%
% exclude control group for analysis
Delta= Delta(default~=0);
donation= donation(default~=0);
default= default(default~=0);

% summary of delta
maxD=max(Delta)
minD=min(Delta)
meanD=mean(Delta)
medianD=median(Delta)

%%
x=[1:300]';
def10=10*ones(length(x),1);
def20=20*ones(length(x),1);
def50=50*ones(length(x),1);

Delta10=Delta_fun(x,def10);
Delta20=Delta_fun(x,def20);
Delta50=Delta_fun(x,def50);

lambda2=lambdas(2);
lambda3=lambdas(3);
mean_delta=(1/lambda2)*ones(length(x),1);
median_delta=log(2)*(1/lambda2)*ones(length(x),1);
%% intersections
Del10_0=@(x1) (x1^2+10^2)/2 - x1*10 - log(2)*(1/lambda2);
Del20_0=@(x2) (x2^2+20^2)/2 - x2*20 - log(2)*(1/lambda2);
Del50_0=@(x3) (x3^2+50^2)/2 - x3*50 - log(2)*(1/lambda2);

zero10=fzero(Del10_0, 190);
zero20=fzero(Del20_0, 200);
zero50=fzero(Del50_0, 210);
Del10=@(x1) (x1^2+10^2)/2 - x1*10
Del20=@(x2) (x2^2+20^2)/2 - x2*20
Del50=@(x3) (x3^2+50^2)/2 - x3*50

%%
figure
plot(x, Delta10, 'k:', 'linewidth',1.2 );
hold on
plot( x, Delta20,'k--','linewidth', 1.2);
plot(x, Delta50, 'k-','linewidth', 1.3);
plot(x, mean_delta, 'k-')
plot(x, median_delta, 'k-.' )
plot(zero10, Del10(zero10), 'ko');
plot(zero20,Del20(zero20), 'k*');
plot(zero50, Del50(zero50),'s', 'MarkerSize',5, 'MarkerEdgeColor','k');
xlabel('\rho','fontweight','bold','fontsize',22) % x-axis label
ylabel('Utils','fontweight','bold','fontsize',18) % y-axis label
legend({'\Delta(\rho,d=10)','\Delta(\rho,d=20)','\Delta(\rho,d=50)', 'mean \delta','median \delta',strcat('rho=',num2str(zero10(1))), strcat('rho=',num2str(zero20(1))), strcat('rho=',num2str(zero50(1)))}, 'Location','northwest','fontsize',12)

h=gcf;
set(h,'PaperOrientation','landscape');
set(h,'PaperUnits','normalized');
set(h,'PaperPosition', [0 0 1 1]);
print(gcf, '-dpdf', 'Delta.pdf');
